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Abstract. A detailed analysis of the century long visual light curve of the long-period Mira star R Cygni is 
presented and discussed. The data were collected from the publicly available databases of the AFOEV, the 
BAAVSS and the VSOLJ. The full light curve consists of 26655 individual points obtained between 1901 and 
2001. The light curve and its periodicity were analysed with help of the O — C diagram, Fourier analysis and 
time- frequency analysis. The results demonstrate the limitations of these linear methods. The next step was to 
investigate the possible presence of low-dimensional chaos in the light curve. For this, a smoothed and noise-filtered 
signal was created from the averaged data and with help of time delay embedding, we have tried to reconstruct 
the attractor of the system. The main result is that R Cygni shows such period-doubling events that can be 
interpreted as caused by a repetitive bifurcation of the chaotic attractor between a period 2T orbit and chaos. 
The switch between these two states occurs in a certain compact region of the phase space, where the light curve 
is characterized by ~1500-days long transients. The Lyapunov spectrum was computed for various embedding 
parameters confirming the chaotic attractor, although the exponents suffer from quite high uncertainty because of 
the applied approximation. Finally, the light curve is compared with a simple one zone model generated by a third- 
order differential equation which exhibits well-expressed period-doubling bifurcation. The strong resemblance is 
another argument for chaotic behaviour. Further studies should address the problem of global flow reconstruction, 
including the determination of the accurate Lyapunov exponents and dimension. 
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1. Introduction 

Red giants on the Asymptotic Giant Branch (AGB) of 
the Hertzsprung-Russell diagram are pulsationally unsta- 
ble and show characteristic light variations. The tradi- 
tional classification (Kholopov et al. 1985-88) is based on 
two properties of the visual light curve, the periodicity 
and full amplitude. Monoperiodic stars with amplitudes 
larger than 2™5 are the Mira-type variables, smaller am- 
plitude or complex light curves place a star among the 
semiregular variables. Even the monoperiodic Mira stars 
do not have strictly repeating light variations, as there 
are apparently irregular cycle-to-cycle changes in ampli- 
tude and/or period. The overwhelming majority of related 
studies addressed the problem of period change, from the 
early works of Eddington & Plakidis (1929), Sterne & 
Campbell (1936) to more recent studies of Isles & Saw 
(1987), Lloyd (1989), Percy & Colivas (1999) and long se- 
ries of papers by Koen & Lombard (from Koen & Lombard 
1993 to Koen & Lombard 2001). Besides a few detections 
of significant gradual period change (e.g. Gal & Szatmary 
1995, Sterken et al. 1999), most of the studies concluded 
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that random period fluctuations dominate the Mira light 
curves. Interestingly, the amplitude variations remained 
quite unstudied, likely because of the nature of the only 
existing data - low precision visual observations. Very few 
investigations can be found on this issue, the most impor- 
tant ones include an analysis of observations by Canizzo 
et al. (1990) and a theoretical study by Icke et al. (1992). 
Both papers inspected the possible presence of chaos in 
long-period variables. While Canizzo et al. (1990) rejected 
the hypothetic low-dimensional chaos in the three stars 
studied, Icke et al. (1992) argued that there is a lot of 
evidence for chaotic behaviour in the red giant variations. 

Variable stars are primary targets of nonlinear analyses 
attempting to detect low-dimensional chaos in astrophys- 
ical systems. The most successful detections are associ- 
ated with pulsating white dwarfs (e.g. Goupil et al. 1988, 
Vauclair et al. 1989), W Vir model pulsations (Buchler 
& Kovacs 1987, Serre et al. 1996b) and two RV Tauri 
type stars, R Set (Kollath et al. 1990, Buchler et al. 
1996) and AC Her (Kollath et al. 1998). A common fea- 
ture is the period-doubling bifurcation, which occurs in 
many simplified model oscillators, too (e.g. Saitou et al. 
1989, Seya et al. 1990, Moskalik & Buchler 1990). Most 



Table 1. A summary of the analysed datasets 
(MJD=JD-2400000). 
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Source 


MJD(start) 


MJD(end) 


No. of points 


AFOEV 


22805 


52181 


8835 


BAAVSS 


15513" 


51907 


11024 


VSOLJ 


22597 


51539 


5888 


VSNET 


49718 


52237 


1105 


Total: 


15513 


52237 


26852 


Corrected: 


15513 


52237 


26655 



a 8 points between MJD 11737 and 13411 were excluded. 

recently, Buchler et al. (2001) presented an overview of 
observational examples of low-dimensional chaos includ- 
ing the first preliminary results for four semiregular stars 
(SX Her, R UMi, RS Cyg and V CVn). The emerging 
view of irregular pulsations is such that large luminosity- 
to-mass ratio strongly enhances the coupling between the 
heat flow and the acoustic oscillation. Consequently, the 
relative growth rates for relevant pulsation modes are of 
order of unity which is the necessary condition for chaotic 
behaviour. Since it is naturally satisfied in the low-mass, 
high-luminosity AGB stars, they are likely candidates for 
chaotic pulsations. 

In this study, we analyse the visual light curve of the 
Mira star R Cygni (= HD 185456 = IRAS 19354+5005, 
spectral type S6/6e with Tc-lines, Jorissen et al. 1998). 
It is a bright, well-observed variable star with an average 
period of ~430 days and average light extrema at ~7? 1 
and ~14™0. The light variation was discovered by Pogson 
in 1852 and the available magnitude estimates date back 
to the late 19th century. It has been studied by numerous 
authors (the ADS lists 179 references), mostly spectro- 
scopically. Wallerstein et al. (1985) pointed out the corre- 
lation between brightness at maximum and interval from 
the previous cycle (fainter maxima occur later than nor- 
mal). R Cyg was also included in a sample of 355 long 
period variables studied by Mennessier et al. (1997) who 
presented a classification of the light curves of LPVs for a 
discrimination between carbon- and oxygen-rich stars. To 
our knowledge, there is no further study dealing with the 
light variation of this star. 

The main aim of this paper is to present a thorough 
analysis of the hundred-years long light curve. The paper 
is organised as follows. The observations and data prepa- 
ration are described in Sect. 2. An intention to interpret 
the light curve with standard linear methods is discussed 
in Sect. 3. The nonlinear analysis is presented in Sect. 4, 
while concluding remarks and further directions are listed 
in Sect. 5. 

2. Observations and data preparation 

Data were taken from the databases of the following or- 
ganizations: the Association Franchise des Observateurs 
d'Etoiles Variables (AFOEVQ), the British Astronomical 




Fig. 1. The distribution of the magnitude differences com- 
puted for 10-day means of the BAAVSS and AFOEV data 
(top) and the BAAVSS and AFOEV data (bottom). The 
maxima at indicate good agreement between the differ- 
ent data. 



Association, Variable Star Section (BAAVSS^) and the 
Variable Star Observers's League in Japan (VSOLjQ). 
Since these data end in early 2001, the latest part of the 
light curve is covered with data obtained via the VSNET 
computer service. 

The basic properties of the datasets are listed in Table 
1. The three main sources are highly independent, as 
very few (less then 1%) duplicated and no triplicated 
data points were found. The total number of estimates 
(duplicates counted twice) is 26852. Before merging the 
data we have checked the data homogeneity with a sim- 
ple statistical test. We have calculated 10-day means for 
all of the three major data sources (BAAVSS, AFOEV, 
VSOLJ) and compared their values in the overlapping re- 
gions. Since the British data are the longest one, we com- 
puted magnitude differences in sense of (BAAVSS) minus 
(AFOEV) and (BAAVSS) minus (VSOLJ). Their distri- 
bution is plotted in Fig. 1. While the AFOEV data result 
in a quite symmetric distribution, there is a suggestion for 
the VSOLJ data for being slightly brighter in average than 
the BAAVSS data. However, the highest peak at Amag=0 
and the symmetric nature for both distributions suggest 
that there is no need to shift any of the data and thus, we 
simply merged all data. After merging the original data 
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tittp : //www.kusastro .kyoto-u. ac . jp/vsnet/gevs 



Kiss & Szatmary: Chaotic behaviour of R Cygni 



3 




10 12 
visual magnitude 

Fig. 2. Top panel: The evolution of number of points per 
bin along the light curve. Middle panel: The same param- 
eter against the apparent magnitude. Note the slightly tri- 
angular distribution with more data in the brighter states. 
Bottom panel: The mean deviation of estimates from the 
mean values (denoted with a). 



files, we checked the deviant points by a close visual in- 
spection of the combined light curve. Almost 200 points 
were rejected and the finally adopted set contains 26655 
visual magnitude estimates. 

The first step in the data processing was averaging the 
light curve using ten-day bins. With the binning procedure 
some statistical findings were revealed. First, we plotted 
the number of points per bin against time (Fig. 2, top 
panel). Here one can find fairly high number of points per 
bin (with steady increase in time) which is a necessary 
condition for calculating accurate mean values. Second, 
the distribution of the same parameter against the appar- 
ent magnitude (Fig. 2, middle panel) shows that there are 
typically 5-10 points per bin even in the faintest state, 
which is enough for good mean values. Finally, the mean 
deviation of the individual points against the apparent 
magnitude (Fig. 2, bottom panel) does not show signifi- 
cant trends, therefore, from a statistical point of view, the 



calculated mean light curve can be regarded quite homo- 
geneous. The typical uncertainty of a single observation 
can be thus estimated as ±0^5. Consequently, the un- 
certainty of the calculated mean values (« ±0™5 / V-^obs) 
ranges from ±0™08 to ±0™17. The binned light curve is 
presented in Fig. 3. The fact, that there are quite long 
subsegments, in which alternating maxima are obviuously 
present, reminded us of the case of RV Tauri type stars, 
i.e. pulsating yellow giants with alternating minima. This 
alternation in R Set and AC Her was interpreted as being 
caused by low-dimensional chaos (Kollath 1990, Kollath 
et al. 1998) and that is why we tackled the question of 
possible chaotic behaviour in the light curve of R Cyg. 
The light curve shape also closely resembles that of some 
pulsating white dwarfs (e.g. Goupil et al. 1988), though 
with a characteristic time-scale of five order of magnitudes 
larger. 

Since we intended to study the light variation with 
tools of the nonlinear time-series analysis, further pre- 
processing was done by smoothing and interpolating the 
binned light curve. First, we fitted Akima splines (Akima, 
1970) to the binned data to get an equidistant series of 
light curve points, keeping the 10-day length (it means 
w43 points per cycle which is a fairly good sampling). The 
spline fit was further smoothed with a Gaussian weight- 
function, where we choose a wider (FWHM=20 days) 
Gaussian than in the case of studies of semiregular vari- 
ables (Kiss et al. 1999, Lebzelter & Kiss 2001). We have 
experimented with various FWHMs starting from 5 days 
up to 30 days; the nonlinear analysis showed that much 
clearer images of the attractor can be reached if the ob- 
servational noise, as well as remaining noise that is ex- 
traneous to the star (e.g. high-frequency components due 
to low- amplitude, but high-dimensional jitter, see Buchler 
et al. 1996) is smoothed out. The adopted 20 days pro- 
vided well-smoothed signal with no zig-zags in the light 
curve and avoided strong amplitude decrease due to over- 
smoothing as in the case of FWHM=30 days. The final 
result is a dataset {s(t n )} that is sampled at constant 10 
days intervals. In Fig. 4 we show some typical observed 
light curve segments with the smoothed and noise-filtered 
signal. Obviously, some sharp features are missed in the 
latter data, however, as will be presented later, the success 
of our nonlinear analysis supports our data processing. In 
order to enable for the reader to study this interesting 
star with other nonlinear methods, we make the original 
10-day means and the smoothed and noise-filtered data 
publicly available in electronic form at CDS, Strasbourg^ 



3. Standard methods 

In this Section the period variation is examined with the 
O — C diagram and one of its descendants, then we show 
that Fourier-analysis fails to give a full description of the 



4 Corresponding data files can be found at CDS via anony- 
mous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via 
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Fig. 3. The whole dataset of R Cygni (10-day means). Note the presence of alternating maxima, especially in the 
lowest two panels. 
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Fig. 4. Typical observed light curve segments (light dots) 
with the smoothed and noise-filtered signal (solid line). 



light curve. Furthermore, fairly strong subharmonic com- 
ponents occur in the frequency spectrum which will be 
naturally included in the adopted nonlinear description 
discussed in the next Section. Two time-frequency meth- 
ods are used to draw some constraints on the stability of 
the frequency content. 



3.1. Period change 

Traditional variable star research studying single-periodic 
light variations strongly relies on the O — C method, al- 
though it has been shown in many cases that Mira stars 
usually do not fit the requirements of its application. This 
is due to the large intrinsic period fluctuations found 
in these variables (see, e.g. recent studies by Sterken et 
al. 1999, Percy & Colivas 1999), but even white noise 
(Koen 1992) or low-dimensional chaotic systems (Buchler 
& Kollath 2001) can generate curved O — C diagrams lead- 
ing to spurious detections of long-term period changes. 
Nevertheless, despite the many possible interpretations, it 
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Fig. 5. The O - C diagram of R Cyg (£ =MJD 18531, 
P=428 d). The typical uncertainty of a given point is 
about ±10 days. 



might be interesting to show the results of this traditional 
method. 

In order to construct the O — C diagram we have 
determined all times of maxima from the binned light 
curve shown in Fig. 3. This was done by fitting low- 
order (3-5) polynomials to the selected parts of the light 
curves around maxima. Typically we selected ±60-70 
days around a given maximum (this means ~13— 15 points 
of the light curve) and fitted polynomials with different 
orders chosen to give optimal fit (as judged by a visual in- 
spection of the data). The fitted interval slightly changed 
with time as several cycles were covered only partially. 
Thanks to the extraordinarily good global light curve cov- 
erage, no cycle was lost between 1901 and 2001, thus we 
could determine 86 epochs of maximum. A mean period of 
428 days given by the Fourier analysis of the whole dataset 
(next subsection) and the eighth epoch (MJD 18531) was 
used to construct the O — C diagram plotted in Fig. 5. 

Three different "linear" segments appear in the O — C 
diagram which could be explained by two period jumps 
at MJD 23500 and MJD 47200. However, little can be 
said about the astrophysical meaning of these jumps. 
Following the work by Percy & Colivas (1999), we per- 
formed the Eddington-Plakidis test (Eddington & Plakidis 
1929) to estimate the amount of random period fluctua- 
tion e. The resulting relative fluctuation e/P is 0.01 that 
places R Cygni close to the lower boundary of Fig. 2 in 
Percy & Colivas (1999), which means a relative stability 
of the mean period compared to other Mira stars of sim- 
ilar periods. Our most important conclusion here is that 
we can exclude the star to be undergoing noticeable long- 
term evolution (due to, e.g., a thermal pulse deep in the 
stellar interior, Wood & Zarro 1981) over the one century 
of observations. The apparent period is shorter at the very 
end of the dataset, but its variation does not seem to be 
a result of a slow period evolution. We note, that Koen & 




o.oi 



0.002 



0.004 0.006 
freq. (c/d) 



0.008 



0.01 



Fig. 6. The frequency spectrum of R Cygni (top panel). 
Note the presence of subharmonic components up to 5/2 
fo- The window function shows two alias peaks, one is the 
1-year alias due to seasonal gaps, while the other one at 
0.00233 d _1 (i.e. 1/430 d _1 ) is caused by a few unobserved 
minima in the light curve. 



Lombard (2001) used simultaneously the epochs of max- 
ima and minima to test for period changes in long-period 
variables which is a more effective approach than that of 
the O — C diagram - at least in those cases, where only 
times of maxima and minima are available. In our case, 
the use of the full light curve is strongly favoured instead 
of examining any particular part of it. 

3.2. Fourier analysis 

To investigate possible multiperiodicity, we have carried 
out a Fourier analysis with subsequent prcwhitening steps. 
We have used Period98 of Sperl (1998) which also includes 
multifrequency least squares fitting of the parameters. The 
studied frequency range was 0-0.01 d" 1 with a step size 
of 1.8 • 10~ 6 . We have also determined the main period in 
all 6 subsets plotted in Fig. 3 to compare apparent period 
change with the suggestions of the O — C diagram. 

The calculated amplitude spectrum of the whole light 
curve is presented in Fig. 6, where the window function 
is shown, too. The primary peak is at /o=0-002338 d _1 
(P = 428 d). The structure of the frequency spectrum 
is quite complex. One can find the integer harmonics of 
the main frequency (2/o, 3/o, 4/o), while subharmonic 
components are present, too (3/o/2, 5/o/2). Note, that 
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Fig. 7. A comparison of the 15 frequency fit with the in- 
terpolated (top) and extrapolated (bottom) data (see text 
for details). The extrapolation obviously fails. 



these subharmonic components are not introduced by the 
Fourier approach through subsequent prewhitening with 
sine curves, as Fig. 6 shows the frequency spectrum before 
any prewhitening. Their presence is due to the alternat- 
ing amplitude changes; this has been proven by a very 
simple test. We have created artificial light curves by re- 
peating i) a single cycle 50 times and ii) a double cycle 
with successive low and high amplitudes 25 times. Then 
the computed Fourier spectra clearly marked the reason 
for the subharmonics. 

The main problem is that instead of well-determined 
narrow peaks, there are closely separated groups of 
peaks at the mentioned frequency values. The first 25 
prewhitening steps resulted in 8 components in the range 
0.002249-0.002464 d^ 1 (« / ), 4 components in the range 
0.004605-0.004781 d" 1 (« 2/ ), 6 components in the 
range 0.001007-0.001437 d" 1 (w / /2), 1 component at 
0.003558 (» 3/ /2), 1 component at 0.007001 (« 3/ ) and 
5 low-frequency components corresponding to the appar- 
ently irregular mean brightness variation. However, even 
the 25-frcqucncy fit leaves a residual scatter of 0™52 being 
much larger than the expected uncertainty of the points. 
Therefore, the fit, besides being physically irrelevant, is 
still not perfect. 

Another test of the multiperiodic hypothesis was done 
as follows. We have taken the first half of data, between 
MJD 15000 and 34000. This subset was Fourier analysed 
with the same procedure. The first 15 frequencies yielded 
an acceptable fit of the whole subset. Then this set of fre- 
quencies (and amplitudes and phases) was used to predict 
the second half of data. As expected, the extrapolation 
failed. This is demonstrated in Fig. 7, where subsets of 
the halves are plotted with the interpolating and extrap- 



olating fits. The extrapolated signal has obviously lost its 
predictive power. Again, we have to reject the multiperi- 
odic description of the light curve. 

The frequency grouping found in R Cyg is exactly the 
same what as the one found for R Set by Kollath (1990) 
and AC Her by Kollath et al. (1998). In those cases de- 
tailed tests showed that multiple periodicity can be ex- 
cluded as the reason for the high number of components 
in the Fourier spectrum. Furthermore, such behaviour is 
characteristic for chaotic systems, e.g., the well-studied 
Rossler oscillator (see Serre et al. 1996a and Buchler 
& Kollath 2001 for detailed discussions in astronomical 
terms). 

The presence of subharmonic components is another 
interesting property. Similar subharmonics were found in 
a number of pulsating white dwarf stars in which they 
were interpreted as indications of period-doubling bifur- 
cation (e.g. Vauclair et al. 1989). The /o/2 subharmonic is 
usually accepted indicator of period-doublings. This sub- 
harmonic is fairly strong in the Rossler oscillator (see, e.g., 
Fig. 3 in Serre et al. 1996a) or other chaotic systems where 
pronounced period-doubling occurs (a closely related ex- 
ample is the case of W Vir hydrodynamic models analysed 
by Serre et al. 1996b). 

The separate analysis of 6 segments, each being ^6000 
days long illustrates the essentially insignificant changes 
of the main period. The resulting periods are the follow- 
ing: MJD 15500-21600: 423±1.5 d; MJD 21600-27700: 
429±2 d; MJD 27700-33900: 428±1 d; MJD 34000-40000: 
428±0.5 d; MJD 40000-46100: 427±3.2 d; MJD 46100- 
52200: 422±1 d. This gives the same picture of the period 
change as the O — C diagram does: slightly shorter periods 
at the beginning and end of the data, between them the 
mean period was constant. However, the given differences 
are consistent with the random fluctuations of the period 
at 1% level. 



3.3. Time-frequency analysis 

Modulations and sudden changes of the light curve (i.e. 
time-dependent frequency content) can be studied with 
several time-frequency methods. In variable star research, 
one of the widely used methods is the wavelet analysis 
(Szatmary et al. 1994 and references therein). A more ef- 
ficient method is the use of of Choi- Williams distribution 
(CWD, Choi & Williams 1989) which enhances the sharp 
features in the time-frequency-amplitude space. 

We present both the wavelet map and CWD in Fig. 8. 
They were calculated with the software package TIFRAN 
(Time FRequency ANalysis) developed by Z. Kollath and 
Z. Csubry at Konkoly Observatory, Budapest. As ex- 
pected, the main ridge at / w 0.0025 d _1 is the most 
obvious feature, while its harmonics are visible up to the 
third order. Subharmonic components can be identified in 
the CWD and we also noticed a slight tilt of the ridges 
(most emphasized for 2/o) in the last third of the data. 
The suggested period shortening around MJD 40000 (note 
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Fig. 8. The wavelet map (left) and Choi- Williams distribution (right) for R Cyg. 



the —15500 shift of the horizontal axes in Fig. 8) is in 
good agreement with the results of the Fourier analysis. 
Unfortunately, little can be inferred about the physical 
system behind the light variation. Time-frequency anal- 
ysis, similarly to Fourier analysis, is basically an inter- 
polation method. There are strict limitations and real un- 
derstanding needs a completely different approach. This is 
provided by the methods of nonlinear time-series analysis. 

4. Nonlinear analysis 

During the last two decades, there has been a growing 
interest in nonlinear studies of variable stars. Several ex- 
amples of seemingly irregular light curves were interpreted 
as results of low-dimensional chaotic behaviour. Buchlcr 
and his co-workers published some fundamental papers 
dealing with either theoretical interpretations or analy- 
ses of observations. There has been a number of excellent 
papers and reviews in the astronomical literature on this 
topic (e.g. Perdang 1985, Serre et al. 1996a, Buchler et al. 
1996, Buchler & Kollath 2001), where the reader can find 
definitions and descriptions of the analysing methods. 

Hereafter we assume that the light curve is generated 
by a deterministic dynamical sytem of a physical dimen- 
sion d. In this section, we employ various techniques to 
unfold the multidimensional structure of this system in 
reconstructed phase spaces. The most important method 
is the time delay embedding. We discuss some proper- 
ties of the embedding and the reconstructed attractor of 
the system, based on the smoothed and noise-filtered light 
curve signal {s(t n )}. For this purpose, we have extensively 
used the TISEAN package of Hegger et al. (1999), which 



is a publicly available software package consisting of prac- 
tical implementations of nonlinear time-series methods^} 
In the following we will refer to some routines of this pack- 
age (e.g. mutual, f alse_aearest, svd etc.). Their use and 
basics are very well documented at the cited website. 



4.1. Embedding parameters 

The first issue is the optimal choice of embedding param- 
eters. As has been discussed in, e.g., Buchler et al. (1996) 
different time delays (A) and/or embedding dimension 
(d c ) may result in completely different phase space recon- 
structions. In the case of R Set and AC Her, a small range 
of fraction of the formal period (5-20%) yielded "optimal" 
results. We applied the routine mutual which calculates 
the time delayed mutual information suggested by Fraser 
& Swinney (1986) as a tool to determine a reasonable de- 
lay. The first minimum of the mutual information occured 
at A=13 (note, that henceforth we normalize A with the 
10-days long sampling rate) and a visual inspection of de- 
lay representations supported the optimal A to be between 
10 and 15. The minimal embedding dimension was esti- 
mated by the false nearest neighbour method proposed by 
Kennel et al. (1992). False neighbours in the phase space 
occur when the reconstruction is done with a d e < d and 
points well-separated in higher dimensions get closer due 
to projections into a lower dimension space. The idea is 
to find that dimension in which the embedding results in 
marginal amount of false neighbours, i.e. a "good" dis- 
tribution of points is reached in the reconstructed phase 



See also frittp : / / www . mpipks-dresden . mpg . de/~t isean 
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space. The routine f alse_nearest calculates the fraction 
of false neighbours while changing d e for a given distance 
threshold. The results of this method should be handled 
carefully, because in case of short datasets, it may under- 
estimate the optimal dimension. For short sets and higher 
dimensions one gets few neighbours simply because the 
points sparsely occupy the embedding space. In our case 
d c ~3 or 4 is suggested by this method, as the fraction of 
false neighbours decreases drastically for larger embedding 
dimensions (even for 5 and 6). Therefore, we adopted 10 
to 15 as the reasonable range for A and 3, 4 for d e . 

4.2. Recurrence plots 

To identify time resolved structures of the dataset, we 
have determined recurrence plots with various parame- 
ters (with the routine recurr). It is a very simple tool: 
it scans the time series and marks each pair of time in- 
dices whose corresponding delay vectors has distance 
< e for a given e. Therefore, a black dot in the (i, j)-plane 
means closeness. If a periodic signal is examined, then the 
points tend to occur both on the diagonal and parallel 
lines separated from the main diagonal by integer multi- 
plets of the period. These plots can be very well used to 
detect and confine transients in a dataset. Fortunately, the 
recurrence plot is not particularly sensitive to the choice 
of embedding parameters and we found only marginal, 
almost invisible differences when choosing embedding di- 
mension and time delay within reasonable ranges (3-4 for 
d e and 10-15 for A). We show the recurrence plot of our 
data for d c — 3 and A = 13 in (top panel in Fig. 9), whose 
values were chosen without any special consideration. In 
this Figure we have converted (i, j) indices to MJD. 

Several interesting conclusions can be drawn even 
from this very simple plot. A few well determined tran- 
sients occured between MJD 16300-17300, 30700-32500 
and 42700-46300 (with a certain ambiguity of the bound- 
aries). The first one is probably caused by the lack of 
observations and thus uncertain spline interpolation, but 
the other two are real as they are well covered by data. 
Thus one has to be careful when analysing subsets that 
include the transient regions. The periodicity is obvi- 
ous from the dense paralcl lines in the plot. One can 
find large islands where the distance between the paral- 
lel lines jumps to the doubled period. To emphasize such 
details, we plot two subrcgions in the middle and bot- 
tom panels of Fig. 9. The latter one clearly corresponds 
to period-doubled region. By a close inspection of the top 
panel, one can find similarly period-doubled regions, e.g., 
around (20000,23000), (23000,41000) or (41000,41000) in 
<~ 1500 x 1500 wide blocks. Furthermore, it is interesting 
that both two large transients were similar as suggested 
by the clustered points around (30000,45000). It follows 
that the transients must have occured around a particu- 
lar small region in the phase space, therefore, they were 
caused by the same process. 
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Fig. 9. The recurrence plot and two subregions with dif- 
ferent behaviour (marked with thick boxes in top panel). 
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Fig. 10. BK projections (d c = 4, A = 10). See text for further detais. 



4.3. Visualization of the trajectories 

To visualize the embedded state vectors s n , we have gener- 
ated the Broomhcad-King (BK) projections (Broomhead 
& King 1986), which project onto the eigenvectors of the 
correlation matrix. In the TISEAN package, this can be 
done with the routine svd (i.e. singular value decomposi- 
tion). We plot the most informative BK projections ob- 
tained from different data. The first column (R Cyg (all)) 
is based on the whole {s(t n )}, with d c = 4, A = 10. The 
second column (R Cyg (40-52)), with the same embedding 
parameters, contains the reconstructed attractor from the 
last third of the data, between MJD 40000 and 52200. 
Plots in the third column (syn (40-52)) are based on a 
synthetic dataset, which has been calculated by a locally 
linear prediction using the routine nstep. The determined 
local linear model can be iterated for an arbitrarily long 
time thus allowing the creation of synthetic datasets. In 
our case, the locally linear approximation was iterated for 
2400 points (i.e. 24000 days). The reason for using only the 
last third of the data for the local linear model is the fact, 
that this way we only approximate the unknown map f of 
the system (Buchler & Kollath 2001) with its local Taylor 
expansion. For longer data the linear model oversmooths 
the trajectories in the phase space thus blurring out the 
phase portrait of the attractor. The fourth column in Fig. 
10 shows plots based on a simple one zone model discussed 
in the next section. 



Our main conclusions are as follows. The embedding 
produces trajectories with remarkable structures. The 
global structure does not depend whether the whole set 
or some of its subsets are analysed. By iterating a locally 
linear approximation of the map we could get a cleaner 
image of the attractor which shows the period-doubling 
unambiguously. It is suggested that the system switches 
back and forth between a period 2T orbit and chaotic 
state. The switch occurs in a located region of the phase 
state and there shows the light curve its transients. 

4.4. Lyapunov exponents and dimension, the 
correlation integral 

The exponential growth of infinitesimal perturbations 
from which the chaos arises can be quantitatively char- 
acterized by the spectrum of Lyapunov exponents (e.g. 
Kantz & Schrciber, 1997 and references therein). If chaos 
is present in a dataset then one of the Lyapunov exponents 
should be positive. We estimated this maximal exponent 
by computing 

S(e,d e ,t) = ( In I pji-r J2 |s n +t - Sn'+t| ) 
\ V s -' elx " //n 

with the routine lyap_k (U n is the e-neighbourhood of s n , 
excluding s n itself). If S(e,d e ,t) shows a linear increase 
with the same slope for all d c larger than some d and for 
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Fig. 11. Estimating the maximal Lyapunov exponent. 
The solid line is the linear fit shifted vertically for clar- 
ity. The slope is 0.0244 implying Ai=0.00244 d" 1 . 



a range of e, then this slope is the estimate of the maxi- 
mal exponent Ai. We calculeted S(e, d e , t) using the whole 
{s(t n )} and plot the result in Fig. 11. Different branches 
corresponding to e? c =3,4,5 and 6 and three values of e 
show oscillations but are generally parallel. A global lin- 
ear fit gives a slope 0.0244 corresponding to Ai =0.00244 
d _1 . This roughy equals to that of determined for R Set 
(«0.0020 d _1 , Buchler et al. 1996) and smaller by a factor 
of 2-3 as for AC Her (Kollath et al. 1998). 

The full Lyapunov spectrum was computed with the 
routine lyap_spec. Its results should be considered as only 
preliminary, because the method, similarly to nstep, em- 
ploys local linear fits. It fits local Jacobians of the lin- 
earized dynamics and by multiplying them one by one to 
different vectors in tangent space, emerges the trajecto- 
ries. The logarithms of rescaling factors are accumulated 
and their average give the Lyapunov exponents (Sano & 
Sawada 1985). 

In order to give an indication of the robustness of the 
Lyapunov spectrum and Lyapunov dimension, we present 
their values for various embedding parameters in Table 2. 
In every cases, we find one positive exponent. Thus, the 
attractor is clearly chaotic. In most cases, the absolute 
value of the second exponent is smaller than the first by 
a factor of 3 to 5. An ideal flow should have A2 = (Serre 
et al. 1996a), therefore our local fit is indeed only a rough 
approximation. Consequently, the determined Lyapunov 
dimensions scatter between 2 and 3 with no definite ten- 
dency. Unfortunately, as noted, e.g. in Serre et al. (1996a), 
the computation of all of the Lyapunov exponents requires 
very long datasets. This is the other reason why our results 
are fairly uncertain. A global nonlinear fit would yield re- 
sults that are superior to ours presented here. We have 
to conclude that our methods do not allow an accurate 
determination and a further study should be devoted to 
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i) a global flow reconstruction (Serre et al. 1996a) and ii) 
computation of more accurate Lyapunov exponents and 
dimension. The latter one is important because these num- 
bers are the main quantitative parameters of a chaotic 
system. 

Finally, the correlation dimension was estimated by the 
correlation sum C(d e ,e) (Grassbergcr & Procaccia 1983). 
This method was applied by Canizzo et al. (1990) for three 
long-period variable stars. In those cases, Canizzo et al. 
concluded that the light curves can be modelled with a 
regular and a stochastic component thus rejecting the hy- 
pothetic chaos. Buchler et al. (1996) argued against the 
use of correlation integral, claiming that typically avail- 
able data are too short for this purpose. Nevertheless, we 
have computed the correlation sum to check its properties. 
For this, the routine c2naive was used which output was 
processed by the routine c2d. The latter routine computes 
dlogC(e)/dloge for which one expects a wide constant 
local minimum with different embedding dimensions (see 
detailed tests presented in Canizzo et al. 1990). 

We show the results of our computations in Fig. 12. 
This diagram is based on the whole {s(t n )}. Although 
we could not find the expected linear scaling region, the 
presented behaviour slightly resembles that of the Lorenz- 
attractor as shown in Canizzo et al. (1990) (see especially 
their Fig. 4b). The closely separated local minima suggest 
a fractal dimension between 2.0 and 2.3 for the attractor, 
though this estimate should be taken with caution. We 
conclude that the analysis of the correlation sum yields 
results that are compatible with the previous ones. 

4.5. Comparison with a one zone model 

We compare the presented behaviour of R Cygni with a 
simple one zone model with stressed period-doubling bi- 
furcation. The similarity is provided by the reconstructed 
attractor. Although it has only marginal physical implica- 
tions, the resemblance is quite illustrative. 
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Fig. 12. Correlation dimension estimation from the cor- 
relation sum. Embedding dimensions 3 to 25 are shown 
for odd numbers. Slopes are determined by straight line 
fits to the log-log plot of the correlation sum. Note that 
horizontal axis is logarithmic. 



As a simple model for the dynamic of a pulsating white 
dwarf star, Goupil et al. (1988) used the following third 
order differential equation 

x'" + Kx" + x' + K(j,x(l + fix) = (1) 

corresponding to a one zone model that is analogous to 
Baker's one zone model (Baker 1966). Here x is the radial 
displacement from equilibrium. The control parameter is 
fi, while K — 0.5 and (3 — —0.5 was fixed (see Goupil et 
al. 1988 for their meaning). When \i is varied, the system 
first bifurcates from its stable fixed point to a stable limit 
cycle (period T orbit at [i — —1), while at [i — — 1.66 the 
first period-doubling bifurcation occurs. The period T or- 
bit becomes mildly unstable and a stable period 2T orbit 
exists. This was the main point which reminded us to the 
case of R Cygni and that is why we have solved Eq. (1) nu- 
merically During the fourth order Rungc-Kutta numeri- 
cal integration we have added simple random perturbation 
that mimicked internal perturbations of the system. To al- 
low an easy comparison with R Cygni, we have rescaled 
the time to get 430 days for the period of the one zone 
model. 

The solution is compared with the light curve and its 
iterated locally linear approximation in Fig. 13. The am- 
plitude does not change as much as in the original data, 
which is simply due to the fact, that Eq. (1) gives the 
variations of the radius, not the luminosity. There is, of 
course, some functional dependence between them, which 
is not expected to change the very basic nature of the 
reconstructed attractor. The calculated maximum alter- 
nation, its disappearance and reappearance is similar to 
what characterizes R Cygni, as has already been stressed 
by the BK projections. A further support is given by the 




MJD (-40000) !3000 

Fig. 13. A comparison of R Cyg (40-52), syn (40-52) sets 
with the one zone model. 



Fourier spectra of data shown in Fig. 14. As discussed in 
Goupil et al. (1988), depending on how long and at what 
time we observe the system, we see an oscillation whose 
trajectory in phase space lies either on the period 2T or- 
bit or on the period T orbit (then disappears the /o/2 
subharmonic), or a mixed trajectory is observed switch- 
ing over from one attractor to the other. In this case, the 
temporal behaviour exhibits a change in the alternation 
of small and large maxima. This largely reminds us of the 
light curve behaviour of R Cyg. 

5. Conclusions 

Our study concentrated on the demonstration of the low- 
dimensional chaotic behaviour in the light curve of a Mira 
star. For this, various linear and nonlinear methods has 
been applied. While the O — C diagram is compatible 
with either stochastic (high dimensional random process) 
or chaotic (low-dimensional deterministic) dynamics, the 
frequency and time-frequency analyses raised the possi- 
bility of underlying chaos through the presence of subhar- 
monic components. The strongest pieces of evidence came 
from the nonlinear analysis. Time delay embedding and 
different analyses of embedded state vectors in the recon- 
structed phase space revealed the following results: 

1. The phase portraits show remarkably regular struc- 
tures. The regularity does not depend whether the full 
dataset or some subsets are used for embedding. The esti- 
mated minimum embedding dimension was 3 or 4 as sug- 
gested the false nearest neighbours test. Other extensively 
used numerical methods favoured 4. 

2. The long term behaviour is dominated with apparent 
period-doubling events. It can be interpreted as caused by 
a chaotic system which switches back and forth between a 
period 2T orbit and chaotic state. The switch occurs in a 
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Fig. 14. Fourier spectra of the data plotted in Fig. 13. 



compact region of the phase space in which the light curve 
exhibits 1000-1500-days long transients. 

3. The Broomhead-King projections were used to visualize 
the phase portraits. Synthetic data were generated with a 
locally linear fit to the attractor. They gave the clear- 
est phase portraits showing unambiguously the period- 
doubling bifurcation. 

4. Quantitative parameters, such as the Lyapunov expo- 
nents and dimension, were estimated, though with fairly 
large uncertainty. The estimated largest Lyapunov expo- 
nent is +0.0024 d _1 . This means that the distances be- 
tween the neighbouring trajectories increase by e during 
417 days, i.e. we cannot predict the variation for longer 
time than one pulsation period. The Lyapunov spectrum 
was computed with the locally linear approximation with 
various embedding parameters. In every cases we found 
one positive exponent. In most case, the second one is close 
to which is characteristic for ideal flows. The correlation 
sum, though rather inconclusively, suggests a fractal di- 
mension 2.0-2.3 for the attractor. 

5. We have compared the analysed data with a simple one 
zone model generated by a third-order differential equa- 
tion. It shows such period-doubling bifurcation which is 
very similar to what is observed for R Cyg. The simi- 
larity of these two systems is our last argument for low- 
dimensional chaos in the light curve of R Cyg. 



How can we compare our results with recent theoret- 
ical advancements in this field of variable star research? 
Most importantly, presence of low-dimensional chaos in a 
Mira star has been detected for the first time. This means 
that despite the complex structure of an AGB star, the 
pulsations in R Cyg occur in a very dimension reduced 
phase space. Presently, we cannot draw firm conclusions 
on the physical dimension of the system but our results 
suggest 3 or 4. In the latter similar consideration 

might be applicable as for R Set by Buchler et al. (1996). 
In that case, they concluded that the pulsations are the 
result of the nonlinear interaction of two vibrationally nor- 
mal modes of the star. As soon as more sophisticated 
nonlinear methods are used to analyse the light curve of 
R Cyg, we will be able to infer more physical implica- 
tions on the nature of pulsations. A further study will 
be devoted to a nonlinear global flow reconstrucion which 
enables, e.g. computing accurate Lyapunov spectrum and 
dimension. Only then it will be possible pinpointing the 
effective (physical) dimension of the phase space. 

The examination of publicly available databases 
(AFOEV, VSOLJ, BAAVSS) revealed that there are at 
least a few dozen or even more than one hundred Mira 
stars with continuous, densely sampled visual light curves. 
R Cygni is only one of them. We list some Mira stars 
of similar period (400-470 days) with well-sampled light 
curves: T Cas, R And, R Lep, R Aur, U Her, \ Cyg, U Cyg, 
V Cyg, S Cep, RZ Peg, R Cas and Y Cas. A quick look 
at their AFOEV and VSOLJ data reveals a few examples 
with interesting light curves resembling that of R Cyg. 
Another further task is to spot those stars which might 
be twins of R Cyg and make an intercomparison between 
them. It is expected that the discovery of other similar 
stars may help to extend the observational base of a new 
discipline, the nonlinear asteroseismology (Buchler ct al. 
1996) which extracts quantitative information from the ir- 
regular light curves. It is, of course, still in its infancy and 
more theoretical as well as observational efforts have to be 
undertaken. 
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